%--- main parameters
% SI System
%clear all;
function [] = Complete_Matlab(option)
global Ux Uy  h  g  p  q  dx dy  Nx Ny  L tmax y;

dx1=1/dx;
dx2 = 1/(dx*dx);

dy1= 1/dy;
dy2 = 1/(dy*dy);


% ---- Construct the spatial Discretization Matrix------%

%------------- Construct the Hyperbolic Term--------------%
A = zeros(Nx*Ny,Nx*Ny);
% Inner Boundary points
for j=2:Nx-1
    for i=2:Ny-1
        index = (j-1)*Ny+i;
        A(index, index-1)= -Uy*0.5*dy1;
        A(index,index+1) = Uy*0.5*dy1;
        A(index, index+Ny) = Ux*0.5*dx1; 
        A(index,index-Ny) = -Ux*0.5*dx1 ;
    end
end

% East and West Boundary
for i=2:Ny-1
    j=1;
    index = (j-1)*Ny+i;
    A(index, index-1)= -Uy*0.5*dy1;
    A(index,index+1) = Uy*0.5*dy1;
    A(index, index+Ny) = Ux*0.5*dx1; 
    A(index,(Nx-1)*Ny+i) = -Ux*0.5*dx1 ;

    j=Nx;
    index = (j-1)*Ny+i;
    A(index, index-1)= -Uy*0.5*dy1;
    A(index,index+1) = Uy*0.5*dy1;
    A(index, i) = Ux*0.5*dx1; 
    A(index,index-Ny) = -Ux*0.5*dx1 ;
end
    
% SW Corner point
i=1;
j=1;
index = (j-1)*Ny+i;
A(index, index+Ny-1)= -Uy*0.5*dy1;
A(index,index+1) = Uy*0.5*dy1;
A(index, index+Ny) = Ux*0.5*dx1; 
A(index,(Nx-1)*Ny+i) = -Ux*0.5*dx1 ;

% SE Corner Point
j=Nx;
index = (j-1)*Ny+i;
A(index, index+Ny-1)= -Uy*0.5*dy1;
A(index,index+1) = Uy*0.5*dy1; 
A(index, i) = Ux*0.5*dx1; 
A(index,index-Ny) = -Ux*0.5*dx1; 



% North and South Boundary
for j=2:Nx-1
    i=1;
    index = (j-1)*Ny+i;
    A(index, index+Ny-1)= -Uy*0.5*dy1;
    A(index,index+1) = Uy*0.5*dy1;
    A(index, index+Ny) = Ux*0.5*dx1; 
    A(index,index-Ny) = -Ux*0.5*dx1 ;
    
    i=Ny;
    index = (j-1)*Ny+i;
    A(index, index-1)= -Uy*0.5*dy1; 
    A(index,index-Ny+1) = Uy*0.5*dy1;
    A(index, index+Ny) = Ux*0.5*dx1; 
    A(index,index-Ny) = -Ux*0.5*dx1 ;    
end

% NW Corner Point
i=Ny;
j=1;
index = (j-1)*Ny+i;
A(index, index-1)= -Uy*0.5*dy1; 
A(index,index-Ny+1) = Uy*0.5*dy1;
A(index, index+Ny) = Ux*0.5*dx1; 
A(index,(Nx-1)*Ny+i) = -Ux*0.5*dx1 ; 

% NE Corner Point

i=Ny;
j=Nx;
index = (j-1)*Ny+i;
A(index, index-1)= -Uy*0.5*dy1; 
A(index,index-Ny+1) = Uy*0.5*dy1;
A(index, i) = Ux*0.5*dx1; 
A(index,index-Ny) = -Ux*0.5*dx1 ;  
%------------- Construction of the Hyperbolic Matrix Finished--------------%

%------------- Compile the Matrices into Global Matrix L--------------%

L=zeros(2*Nx*Ny,2*Nx*Ny);
L(1:Nx*Ny,1:Nx*Ny)=A;
L(Nx*Ny+1:2*Nx*Ny,Nx*Ny+1:2*Nx*Ny)=A;



%------------- Construction of the Elliptic Term --------------%
A = zeros(Nx*Ny,Nx*Ny);

% Populate the diagonal term
for j=1:Nx;
    for i=1:Ny
        index = (j-1)*Ny+i;
        A(index,index) = -2*h*(dx2+dy2);
    end
end

% Inner Boundary points
for j=2:Nx-1
    for i=2:Ny-1
        index = (j-1)*Ny+i;
        A(index, index-1)= h*dy2;
        A(index,index+1) = h*dy2;
        A(index, index+Ny) = h*dx2; 
        A(index,index-Ny) = h*dx2 ;
    end
end



% East and West Boundary
for i=2:Ny-1
    j=1;
    index = (j-1)*Ny+i;
    A(index, index-1)= h*dy2;
    A(index,index+1) = h*dy2;
    A(index, index+Ny) = h*dx2; 
    A(index,(Nx-1)*Ny+i) =h*dx2 ;

    j=Nx;
    index = (j-1)*Ny+i;
    A(index, index-1)= h*dy2;
    A(index,index+1) = h*dy2;
    A(index, i) = h*dx2; 
    A(index,index-Ny) = h*dx2 ;
end
    
% SW Corner point
i=1;
j=1;
index = (j-1)*Ny+i;
A(index, index+Ny-1)= h*dy2;
A(index,index+1) = h*dy2;
A(index, index+Ny) = h*dx2; 
A(index,(Nx-1)*Ny+i) =h*dx2 ;

% SE Corner Point
j=Nx;
index = (j-1)*Ny+i;
A(index, index+Ny-1)= h*dy2;
A(index,index+1) = h*dy2; 
A(index, i) = h*dx2; 
A(index,index-Ny) = h*dx2; 



% North and South Boundary
for j=2:Nx-1
    i=1;
    index = (j-1)*Ny+i;
    A(index, index+Ny-1)= h*dy2;
    A(index,index+1) = h*dy2;
    A(index, index+Ny) = h*dx2; 
    A(index,index-Ny) = h*dx2 ;
    
    i=Ny;
    index = (j-1)*Ny+i;
    A(index, index-1)= h*dy2; 
    A(index,index-Ny+1) = h*dy2;
    A(index, index+Ny) = h*dx2; 
    A(index,index-Ny) = h*dx2;    
end

% NW Corner Point
i=Ny;
j=1;
index = (j-1)*Ny+i;
A(index, index-1)= h*dy2; 
A(index,index-Ny+1) = h*dy2;
A(index, index+Ny) = h*dx2; 
A(index,(Nx-1)*Ny+i) = h*dx2; 

% NE Corner Point

i=Ny;
j=Nx;
index = (j-1)*Ny+i;
A(index, index-1)= h*dy2; 
A(index,index-Ny+1) = h*dy2;
A(index, i) = h*dx2; 
A(index,index-Ny) = h*dx2 ; 


%------------- Compile the Matrices into Global Matrix L--------------%
L(1:Nx*Ny,Nx*Ny+1:2*Nx*Ny)=A;

for k=1:Nx*Ny
    L(Nx*Ny+k,k)=g;
end


%Run the Matlab ODE Solver

y0=[p;q];
dt=1;

tspan = 0:dt:tmax;  % time vector,
options1=odeset('RelTol',1e-6,'Stats','on');
options2=odeset('RelTol',1e-6,'Stats','on');

tic
if option ==1
    disp('ODE 45 Explicit Method Selected');
    [t,y]=ode45(@waves_benchmark,tspan,y0,options1);
else
    disp('ODE 15s Implicit Method Selected');
    [t,y]=ode15s(@waves_benchmark,tspan,y0,options2);
    % [t,y]=ode15s(@upwinding,tspan,y0,options2);
    
end
toc


